% Coreg 
% Chris Whalen
% University of Illinois
% 05/20/04 - Date of Birth
% 01/17/05 - Date of Last Modification

% Coreg - the purpose of this program is to align digitized points (ELP) with
%         the surface of the head from a MRI for a 'best fit'

% Input Files: 1. elp file              4. scalp_zip_file (optional)
%              2. tai file              5. mtg file
%              3. acpc.BRIK             6. hsp file (optional)
%                                        7. roi file (optional)
%              (can also accept src file and det file instead of mtg
%               with change of single command line)

% Output Files: 1. TOL (Talairach Optode Location) file      - aligned talairach transformed ASCII file of optode name/label and x, y, and z coordinates
%               2. AFNI fim file                             - file that can be viewed in AFNI acpc view which shows the location of fiducials and optodes (can also display dummy aligning points by request)
%               3. Segmented Scalp



% When using ELP and HSP to align - the first alignment is done using the
% HSP and subsequent ones are done using ELP. Hence, one should always
% cycle through dbap at least 2x. NB - the second alignment may be worse
% than the first.

clear all
clc

% clc
pause(0.1)
fprintf('\n      COREG    \n\n');

% Initial Parameter Input
data_path=input('Is your data path C:\data\? (y/n) ','s') ;
if data_path=='y'
    data_path='C:\Data\';
else
    data_path=input('Enter your data path: ','s');
end

analysis_path=input('Is your analysis path C:\analysis\matlab\ (for align2 & coreg folders)? (y/n) ','s') ;
if analysis_path=='y'
    analysis_path='C:\analysis\matlab\';
else
    analysis_path=input('Enter your analysis path: ','s');
end

exp=input('Enter experiment folder name: ','s');
prefix=input('Enter filename prefix: ','s');
S=input('Enter subject #: ','s');
date=input('Enter MRI date: ','s');

% Changing Subject Num to Be compatible with dbap
S_dbap = S;
if S(1)=='0'
    S_dbap(1) = '';
    if S(2)=='0'
        S_dbap(1)='';
    end
end

% Dimensions
mri_dim.x = 191; % L -> R
mri_dim.y = 236; % A -> P
mri_dim.z = 171; % I -> S

num_dum_pts = 0; % dummy point initialization
xyz_elp = 0; % initializing
xyz_pre_plot = 0; % initializing
vol = 0;
re = 0; fre = 0; tre = 0;

% Define directories and checking prerequisite files
%
[dir,file] = mkdirs(data_path,prefix,exp,S,date,analysis_path);

choice.mode = input('\nEnter Coreg Mode:\n1 - Standard\n2 - Advanced\n');
check_input_files(file,choice);

% Obtaining Optode Labels and Coordinates
% Align with hsp or elp
choice.align = 2; % Setting Default
if choice.mode == 2
    choice.align = input('\nAlign With: \n1 - Headshape\n2 - ELP\n3 - HSP & ELP\n');
    if (choice.align == 1)
        xyz = hsp2xyz(dir.hsp,exp,[S suffix]);
    end
    if (choice.align == 2)
        [xyz, opt_label] = elp2xyz(file.elp);
    end
    if (choice.align == 3)
        xyz = hsp2xyz(dir.hsp,exp,[S suffix]);
        [xyz_elp, opt_label] = elp2xyz(dir.elp,exp,[S suffix]);
    end
elseif choice.mode == 1
    [xyz,opt_label] = elp2xyz(file.elp); % PRI order; (LPA-RPA)/2 origin
end

% Nasion Asymmetry Correction
if choice.mode == 2
    [xyz, xyz_elp] = fid_asym_correct(choice,xyz,xyz_elp); % Data from 09/07/05 suggests no improvement from applying Nasion Correction
end

% Placing tai fiducials in ref frame w/ tai origin (LPA+RPA)/2 and RAI order
tai = read_tai(file.tai);
tai = tai2brk_ctr(tai,mri_dim);

% Read ROI (Optional)
% roi = []; % initializing
% fid=fopen(file.roi);
% if fid ~= -1 % Determining if ROI file even exists
%     [roi] = read_roi(file,mri_dim); % Reads in ROI file 
%     %fprintf('\nroi.tai_xyz line 31:\n'); % debug
%     %roi.tai_xyz % debug
%     [tai,roi] = tai2acpc(tai,roi); % Under Construction 
%     [roi] = roipts2xyzlist(roi,mri_dim);  % acpc space
%     [roi] = roi2brik(roi,mri_dim);       % making a brik with all faces of roi = 1 and other entries = 0
% %   [roi] = roipts2xyzlist2(roi,mri_dim); % tai space
% %   scatter3(roi.xyz2(:,1),roi.xyz2(:,2),roi.xyz2(:,3),'g','filled'); % debug
% %   hold on;
% %     scatter3(roi.xyz(:,1),roi.xyz(:,2),roi.xyz(:,3),'b','filled');
% %     hold on;
% end

% Placing Polhemus reference system into brk_ctr reference coords (RAI) with mdpt(LPA-RPA) as origin
if choice.align ~= 4 % This extra bit is to check in cap study
	xyz = pol2brk_ctr(xyz);
	if choice.align == 3
        xyz_elp = pol2brk_ctr(xyz_elp);
	end
elseif choice.align == 4
    % permute .tai coords to RAI order
	for i_xyz=1:size(xyz,1)
      temp(1)=mri_dim.x - xyz(i_xyz,2);  % Switching direction of axis 
      temp(2)=xyz(i_xyz,1);  % 1
      temp(3)=xyz(i_xyz,3);  % 3
      xyz(i_xyz,:)=temp;
	end
    xyz = xyz - repmat(tai.origin,length(xyz), 1); 
end

% Generate/Load and View Mask
[choice,file,dir,vol_mask,vol,xyz_mask,thresh] = gen_load_mask_view(choice,file,dir,mri_dim,vol,tai);

% Save Scalp Option - in the future must save as brik format and load as brik file
fid=fopen(file.scalp_zip);
if (fid == -1 | strcmpi(choice.view_mri,'y'))
    save_scalp(file,choice,vol_mask,xyz_mask);
end

mri_markers = [0 0 0]; % initializing
% Aligning Markers to Elp file (CAP study only)
if strcmp(exp,'cap') & choice.mode == 2
    choice.coreg_markers = input('\nWould you like to align marker file to elp file? (y/n)','s');
    if strcmpi(choice.coreg_markers,'y')
        % Doing a sample for dif number of fids
        for num_fids = 32 %3:32
            
            for i=1:length(mri_markers) % Getting 2nd coordinate with correct origin location
                mri_markers(i,2) = mri_dim.x - mri_markers(i,2);
            end
            i=2;j=1;k=3;
            elp_sgn=[1 1 1];
            for i_loc=1:size(mri_markers,1)  % Swap columns 1 with 2
                temp(1)=elp_sgn(1)*mri_markers(i_loc,i);
                temp(2)=elp_sgn(2)*mri_markers(i_loc,j);
                temp(3)=elp_sgn(3)*mri_markers(i_loc,k);
                mri_markers(i_loc,:)=temp;
            end
            %mri_markers = mri_markers - 1;           % Subtracting 1 mm to go back into BRIK coords (?)
            mri_markers = mri_markers - repmat(tai.origin,length(mri_markers),1);
            
            randnums = randperm(32);
            randorder = randnums(1:num_fids); % the random order of points taken to coregister.

            [R,T] = get_lsrbtfm(xyz(3:32,:),mri_markers(3:32,:))
            %[R,T] = get_lsrbtfm(xyz(randorder,:),mri_markers(randorder,:))
            % Transforming all points
            for i=1:size(xyz,1)
                xyz_ls(i,:)=(R*xyz(i,:)' + T)';
            end
            pause;
            dist_mat = sqrt((xyz(:,1)-xyz_ls(:,1))^2 + (xyz(:,2)-xyz_ls(:,2))^2 + (xyz(:,3)-xyz_ls(:,3))^2);
            
            for j=1:num_fids % defining the distance of fiducial point to be zero
                dist_mat(randorder(j)) = 0;
            end
            dist_mat = nonzeros(dist_mat); % removing fiducial points
            
            tre.mean(num_fids) = mean(dist_mat);
            tre.std(num_fids) = std(dist_mat);
            tre.n(num_fids) = length(dist_mat);

            % scatter3(xyz(:,1)-xyz_ls(:,1),xyz(:,2)-xyz_ls(:,2),xyz(:,3)-xyz_ls(:,3) )

%             if strcmpi(choice.view_mri,'y')==1
%                 scatter3(xyz_ls(1:length(mri_markers),1),xyz_ls(1:length(mri_markers),2),xyz_ls(1:length(mri_markers),3),'r','filled');
%                 hold on;
%                 scatter3(mri_markers(:,1),mri_markers(:,2),mri_markers(:,3),'g','filled');
%             end
%             xyz_ls = xyz_ls + repmat(tai.origin,length(xyz_ls),1);
%             [xyz_aol] = mk_aol(xyz_ls,dir,opt_label,exp,S,mri_dim);
        end % for loop num_fid
        plot(tre.mean);
    end % end choice condition
end %



% Coregistering xyz to the tai.fid
if choice.mode == 2 % advanced mode
    fprintf('\nPerform Initial Alignment With:\n1 - Fiducial CW\n2 - Fiducial Least Squares\n3 - Moment Matching\n');
    choice.init_align = input('','s');
    if strcmp(choice.init_align,'1') & choice.align == 3
        fprintf('\nHSP');
        xyz = coreg_fid_cw(xyz,tai);
        fprintf('\nELP');
        xyz_elp = coreg_fid_cw(xyz_elp,tai);
    elseif strcmp(choice.init_align,'1') & choice.align ~= 3
        xyz = coreg_fid_cw(xyz,tai);
    end
    if strcmp(choice.init_align,'2') & choice.align == 3
        fprintf('\nHSP');
        [xyz,fre] = coreg_fid_ls(xyz,tai);
        fprintf('\nELP');
        [xyz_elp,fre] = coreg_fid_ls(xyz_elp,tai);
    elseif strcmp(choice.init_align,'2') & choice.align ~= 3
        [xyz,fre] = coreg_fid_ls(xyz,tai);
    end
    if strcmp(choice.init_align,'3')
        xyz = coreg_moments(xyz,xyz_mask);
    end
elseif choice.mode == 1 % standard mode
    xyz = coreg_fid_cw(xyz,tai); % default is the least-squares method
end
% hold on;
% scatter3(xyz(:,1)+tai.origin(1,1),xyz(:,2)+tai.origin(1,2),xyz(:,3)+tai.origin(1,3),2.5,'r','filled'); % xyz final


% Only Register Non Marker elp points (CAP Study only)
if strcmp(exp,'cap')
    if strcmp(S_dbap,'168')
        test = xyz(18:617,:); 
    else
        test = xyz(33:632,:); % To do just the non-vitE points for cap study: Paradigm 1
    end
else
    test = xyz;
end
ref = xyz_mask;
% % test = xyz(1:32,:);   % to fit just based on vitE points.
% % test = xyz(70:669,:); % To do just the non-Beekley points for cap study: Paradigm 2
% % test = xyz(1:69,:);   % Just Beekley points
% % test = xyz(70:469,:); % No Beekley, no face
% end

% ECC Only: Align Marker File with ELP
if strcmpi(exp,'ecc') 
    choice.ecc = input('\nAlign ECC marker files? (y/n) ','s');
    if strcmp(choice.ecc,'y')
        fprintf('\nECC Study: Aligning MRI IZI Markers\n');
        if ispc
            izi_file = ['C:\data\ecc\reg\' S '_izi.txt' ];
        elseif isunix
            izi_file = ['/data/ecc/reg/' S '_izi.txt' ];
        end
        IZI = importdata(izi_file);
        mri_izi = IZI(:,2:4); % mri_izi coordinates
        mri_izi_tmp = mri_izi; %swapping dimensions
        mri_izi(:,1) = mri_izi_tmp(:,2);
        mri_izi(:,2) = mri_izi_tmp(:,1);
            % Placing into correct reference frame
            for i=1:length(mri_izi) % Getting 2nd coordinate with correct origin location
                mri_izi(i,2) = mri_dim.x - mri_izi(i,2);
            end
            i=2;j=1;k=3;
            elp_sgn=[1 1 1];
            for i_loc=1:size(mri_izi,1)  % Swap columns 1 with 2
                temp(1)=elp_sgn(1)*mri_izi(i_loc,i);
                temp(2)=elp_sgn(2)*mri_izi(i_loc,j);
                temp(3)=elp_sgn(3)*mri_izi(i_loc,k);
                mri_izi(i_loc,:)=temp;
            end
            %mri_markers = mri_markers - 1;           % Subtracting 1 mm to go back into BRIK coords (?)
            mri_izi = mri_izi - repmat(tai.origin,length(mri_izi),1);
        
        
        for i=1:size(IZI,1) % generating elp marker file
            mark_num = IZI(i,1);
            elp_izi(i,:) = xyz(mark_num + 3,:);
        end
        elp_izi
        mri_izi
        [R,T] = get_lsrbtfm(elp_izi,mri_izi)
        %[R,T] = get_lsrbtfm(xyz(randorder,:),mri_markers(randorder,:))
        % Transforming all points
        for i=1:size(xyz,1)
            xyz(i,:)=(R*xyz(i,:)' + T)';
        end
        for i=1:size(elp_izi,1)
            elp_izi_fin(i,:)=(R*elp_izi(i,:)' + T)';
        end
        
        scatter3(xyz(:,1),xyz(:,2),xyz(:,3),'b','filled');
        hold on;
        scatter3(mri_izi(:,1),mri_izi(:,2),mri_izi(:,3),'r','filled');
        scatter3(elp_izi(:,1),elp_izi(:,2),elp_izi(:,3),'g','filled');
        fre_init = mean(sqrt(sum(((elp_izi-mri_izi).^2),2))); % FRE
        fre_fin = mean(sqrt(sum(((elp_izi_fin-mri_izi).^2),2))); % FRE
        fprintf('\nFRE');
        fprintf('\nInitial Average Map Error: %3.1f mm  ', fre_init);
        fprintf('\nFinal   Average Map Error: %3.1f mm\n', fre_fin );
    end
end


% Removing face of back of head points
if  choice.mode == 2
    [test,choice] = rm_face_or_back(test,choice);
else
    choice.rm_pts = 4; % Setting the default to not remove any points
end


% % Randomly Selecting Number of Digitized Points for Registration (CAP study only)
% if strcmpi(exp,'cap') & choice.mode == 2
%     num_test_pts = input(['Enter Number of Test (Polhemus) Points Desired(max ' num2str(size(xyz,1)-32) '):']);
%     num_ref_pts = input(['Enter Number of Ref (MRI) Points Desired(max ' num2str(size(xyz_mask,1)) '):']);
%     test = zeros(num_test_pts,3);
%     ref  = zeros(num_ref_pts ,3);
%     p = randperm(length(xyz)); % giving a random permutation of test points
%     if num_test_pts ~= length(xyz)
%         for i=1:num_test_pts
%             test(i,:) = xyz(p(i),:);
%         end
%     elseif num_test_pts == length(xyz)
%         test = xyz;
%     end
%     q = randperm(length(xyz_mask)); % giving a random permutation of ref points
%     if num_ref_pts ~= length(xyz_mask)
%         for i=1:num_ref_pts
%             ref(i,:) = xyz_mask(q(i),:);
%         end
%     elseif num_ref_pts == length(xyz_mask)
%         ref = xyz_mask;
%     end
% else
%     num_test_pts = 0; % defining to avoid error
% end

    

% Saving 'test' and 'ref' to align2 directory
save_test_ref(dir,S_dbap,test,ref);


% Open Distance Based Alignment Program
%dbap; % If we don't like having 2 Matlab windows we can call dbap from coreg program
                
%scatter3(xyz(:,1),xyz(:,2),xyz(:,3),'b','filled');
%hold on              

% Apply Transforms
num_test_pts = 0;
%[xyz,xyz_elp,choice] = apply_transforms(choice,xyz,xyz_elp,dir,S_dbap,cmp);
[xyz,xyz_elp,choice] = apply_transforms(choice,xyz,xyz_mask,xyz_elp,dir,S_dbap,opt_label,exp,S,mri_dim,mri_markers,num_test_pts,tai);


% Place Digitized Points on Surface of Head
% choice.opt2surf = input('\nForce All Optodes To Lie On Scalp? (y/n) ','s');  
% if strcmpi(choice.opt2surf, 'y') == 1
%     xyz_b4surfforce = round(xyz)
%     [xyz] = opt2surf(xyz,xyz_mask);
%     xyz_aftrsurfforce = round(xyz)
%     xyz_diff = round(xyz_aftrsurfforce - xyz_b4surfforce)
% end

% % Finding Closest Point that is on the surface of the head.
% choice.opt2surf = input('\nForce All Optodes To Lie On Scalp? (y/n) ','s');  
% if strcmpi(choice.opt2surf, 'y') == 1
%     xyz_b4surfforce = xyz;
%     [xyz] = opt2surf2(xyz,xyz_mask);
%     %[xyz] = opt2surf(xyz,xyz_mask);
%     %xyz_aftrsurfforce = round(xyz);
%     xyz_diff = round(xyz - xyz_b4surfforce);
% end



% Place in BRIK coords
xyz = xyz + repmat(tai.origin,length(xyz), 1); 
xyz_mask = xyz_mask + repmat(tai.origin,length(xyz_mask), 1); 



% Obtain src and det from mtg file
% [src,det]=get_src_det(choice.align,src_file,det_file,xyz); % Gets unique src and det and places in arrays
[src,det,mtg]=mtg2src_det(choice,file,xyz); % Gets unique src and det and places in arrays

% % SRC and DET (M. Jacob only)
% if strcmp(exp,'glo') & strcmp(S,'530')
%     fprintf('\nUsing src and det for M. Jacob banana generation\n');
%     src = {'LA01';'LA02';'LA03';'LA04';'LA05';'LA06';'LA07';'LA08';'LB01';'LB02';'LB03';'LB04';'LB05';'LB06';'LB07';'LB08';'LC01';'LC02';'LC03';'LC04';'LC05';'LC06';'LC07';'LC08';'LD01';'LD02';'LD03';'LD04';'LD05';'LD06';'LD07';'LD08'};
%     det = {'LA01';'LA02';'LA03';'LA04';'LA05';'LA06';'LA07';'LA08';'LB01';'LB02';'LB03';'LB04';'LB05';'LB06';'LB07';'LB08';'LC01';'LC02';'LC03';'LC04';'LC05';'LC06';'LC07';'LC08';'LD01';'LD02';'LD03';'LD04';'LD05';'LD06';'LD07';'LD08'};
% %     src = {'LA01';'LA02';'LA04';'LA05';'LA07';'LA08';'LB01';'LB02';'LB04';'LB05';'LB07';'LB08';'LC01';'LC02';'LC04';'LC05';'LC07';'LC08';'LD01';'LD02';'LD04';'LD05';'LD07';'LD08'};
% %     det = {'LA03';'LA06';'LB03';'LB06';'LC03';'LC06';'LD03';'LD06'};
% end

% Inserting Fiducials For HSP
if choice.align == 1   % Labeling the first three sources as fiducials as the fiducials
    opt_label ={'NZ';'LPA';'RPA'};
    for i=4:(length(src)+3)
        opt_label{i} = src{i-3};
    end
end

% Place MRI fiducials in reference frame with (LPA+RPA)/2 as origin
tai.fid = tai.fid + repmat(tai.origin,3,1); 


% Scaling
if choice.mode == 2
    choice.scaling = input('\nScale Data? \n1 - No Scale\n2 - Median Scale\n3 - Regression Scale\n','s');  
    if strcmpi(choice.scaling, '2')
        xyz = scaling(xyz,xyz_mask,mri_dim);
    elseif strcmpi(choice.scaling, '3')
        [xyz,err] = ls_scale(xyz,xyz_mask,mri_dim);
    end
else
    [xyz,err] = ls_scale(xyz,xyz_mask,mri_dim);
end


% Place in BRIK coords
xyz = xyz - repmat(tai.origin,length(xyz), 1); 
xyz_mask = xyz_mask - repmat(tai.origin,length(xyz_mask), 1); 

% Finding Closest Point that is on the surface of the head.
if choice.mode == 2 
    choice.opt2surf = input('\nForce All Optodes To Lie On Scalp? \n1 - No Force\n2 - Closest Point Force\n3 - Radial Projection Force\n');  
    if choice.opt2surf == 2
        % xyz_b4surfforce = xyz;
        [xyz] = opt2surf2(xyz,xyz_mask);
        %[xyz] = opt2surf(xyz,xyz_mask);
        %xyz_aftrsurfforce = round(xyz);
        % xyz_diff = round(xyz - xyz_b4surfforce);
    elseif choice.opt2surf == 3 % radial projection
        [xyz] = opt2surf(xyz,xyz_mask);
    end
elseif choice.mode == 1
    [xyz] = opt2surf2(xyz,xyz_mask);
end

% Creating file to calculate residual error for cap study
if strcmpi(exp,'cap') & choice.mode == 2  
    if strcmp(S_dbap,'168')
        test = xyz(18:617,:); 
    else
        test = xyz(33:632,:); % To do just the non-vitE points for cap study: Paradigm 1
    end
    save ([dir.test 'test' S_dbap '000'], 'test'); % Saving test matrix to go into dbap
    save ([dir.ref  'ref'  S_dbap '000' ], 'ref'); % saving ref matrix to go into dbap
end

% Place in BRIK coords
xyz = xyz + repmat(tai.origin,length(xyz), 1); 
xyz_mask = xyz_mask + repmat(tai.origin,length(xyz_mask), 1); 


% Overplot with Aligned Data
choice.overplot = 'n';
if choice.tfm ~= 1 % Give choice to overplot if transformed the data
    choice.overplot = input('\nOverplot Newly Aligned Optode Positions? (y/n) ','s');  
end
if strcmpi(choice.overplot, 'y') == 1 
    hold on;
    scatter3(xyz(:,1),xyz(:,2),xyz(:,3),1.5,'g','filled');
    rotate3d on;
    % legend('MRI','original (green)', 'aligned (red)');
%     if size(xyz_pre_plot,1) ~= 0
%         for i=1:length(xyz)
%             plot3([(xyz_pre_plot(i,1)) (xyz(i,1))],[(xyz_pre_plot(i,2)) (xyz(i,2))],[(xyz_pre_plot(i,3)) (xyz(i,3))],'k')
%         end
%     end
    %scatter3(xyz(:,1)+tai.origin(1,1), xyz(:,2)+tai.origin(1,2), xyz(:,3)+tai.origin(1,3), 2.5, 'r','filled');
    %scatter3(xyz(1:70,1), xyz(1:70,2), xyz(1:70,3), 2.5, 'r','filled');
end



% Write BRIK
%
% Determining brik suffix
brik_suffix = '';
if choice.align == 1
    brik_suffix = 'hsp';
elseif choice.align == 2
    brik_suffix = 'elp';
elseif choice.align == 3
    brik_suffix = 'elp+hsp';
end
if choice.tfm == 1
    brik_suffix = [brik_suffix '_NoTfm'];
elseif choice.tfm == 2
    brik_suffix = [brik_suffix '_RBTfm'];
elseif choice.tfm == 3
    brik_suffix = [brik_suffix '_Affine'];
end
file.brik_out = [file.brik_out brik_suffix];

% Write acpc vol (BRK if pc, BRIK if unix)
choice.rm_dum = 'n'; % Initializing rm_dum_choice. This is required for later conditions
choice.fbrik = input('\nWrite out optode brik to view in Afni(Create FIM file)? (y/n) ','s');
if strcmpi(choice.fbrik,'y')
    choice.rm_dum = input('Remove dummy ELP points from FIM (if they exist)? (y/n) ','s');
    if ((choice.align == 2) | (choice.align == 3)) & strcmpi(choice.rm_dum,'y') 
        [opt_label,xyz,num_dum_pts]=rm_dum_pts(src,det,opt_label,xyz); % Remove Dummy Points from opt_label and xyz_tai
        if num_dum_pts ~= 0
            file.brik_out = [file.brik_out '_no_dum_pts'];
        end
    end
	fprintf('Display Bananas In FIM:\n1 - No Bananas\n2 - Bananas\n3 - Both (Create 2 FIM files 1 with & 1 without Bananas)\n');
    choice.banana = input(''); % Gives choice to show bananas
	if  (choice.banana == 1) | (choice.banana == 3)  % No Banana FIM
        % Create and Write Optode BRIK
        dig_brik=zeros(mri_dim.x,mri_dim.y,mri_dim.z);  % RAI
        dig_brik = afni_dig_brik(choice,xyz,dig_brik,src,det,tai.fid,opt_label,tai.origin,mri_dim); % Creates the digitized brik with the appropriate colors and parameters
        write_brik(dig_brik,file,dir,exp,S,date);
    elseif (choice.banana == 2) | (choice.banana == 3) % Banana FIM
        file.brik_out = [file.brik_out '_Bananas'];
        fprintf('\nRemoving Dummy Points For Banana Creation ...\n');
        [opt_label,xyz,num_dum_pts]=rm_dum_pts(src,det,opt_label,xyz);
        % Head Mask
        choice.head_mask = input('\nUse Head Mask?\n','s');
        if strcmpi(choice.head_mask,'y')
            [head_mask] = brk2headmsk(file,mri_dim,vol,thresh); % later can make into saved and loadable file
            %view_acpc(head_mask,thresh);
        else
            head_mask = ones(mri_dim.x,mri_dim.y,mri_dim.z); %initiating
        end
        % Brain Mask
        choice.brain_mask = input('\nUse Brain Mask?\n','s');
        if strcmpi(choice.brain_mask,'y')
            [brain_mask] = brk2brainmsk(choice,file,dir,mri_dim)
        else
            brain_mask = ones(mri_dim.x,mri_dim.y,mri_dim.z); %initiating
        end
        [dig_brik] = bananas(choice,file,xyz,src,det,tai,opt_label,mri_dim,head_mask,brain_mask);
        choice.roi_mask = input('\nUse ROI Mask?\n','s');
        if strcmpi(choice.roi_mask,'y')
            dig_brik = dig_brik.*roi.msk; % masking all non-ROI points
        end
        % Banana ROI stats
        fclose(fid);
        fid = fopen(file.roi);
        choice.ban_stats=input('Compute Banana Statistics? (y/n) ','s');
        if (fid~= -1) & strcmpi(choice.ban_stats,'y')
            mri.ban_sum = sum(sum(sum(dig_brik(:,:,:)))); 
            mri.vol = mri_dim.x*mri_dim.y*mri_dim.z;
            roi_brik = cell(roi.num,1); % must define as cells because the roi briks are different size (i.e. can't add to structure).
            roi_brik_reshape = cell(roi.num,1); 
            for i=1:2:2*roi.num
                roi_brik{(i+1)/2, 1} = dig_brik(roi.acpc_xyz(i,1):roi.acpc_xyz(i+1,1),roi.acpc_xyz(i,2):roi.acpc_xyz(i+1,2),roi.acpc_xyz(i,3):roi.acpc_xyz(i+1,3));
                roi.ban_sum((i+1)/2) = sum(sum(sum(roi_brik{(i+1)/2})));
                % summing the banana elements of roi
                roi_brik_reshape{(i+1)/2,1} = reshape(roi_brik{(i+1)/2},[1 size(roi_brik{(i+1)/2},1)*size(roi_brik{(i+1)/2},2)*size(roi_brik{(i+1)/2},3)]); % Converting to single Nx1 vector
                roi.ban_sd((i+1)/2) = std(roi_brik_reshape{(i+1)/2});
                roi.vol((i+1)/2) = (roi.acpc_xyz(i+1,1)-roi.acpc_xyz(i,1))*(roi.acpc_xyz(i+1,2)-roi.acpc_xyz(i,2))*(roi.acpc_xyz(i+1,3)-roi.acpc_xyz(i,3)); % probably contains errors
                roi.ban_mean((i+1)/2) = roi.ban_sum((i+1)/2)/roi.vol((i+1)/2);
                fprintf('\nROI Number: %2.0f', (i+1)/2);
                fprintf('\nROI Banana Sum: %3.1f',roi.ban_sum((i+1)/2));
                fprintf('\nROI Banana Mean: %3.1f',roi.ban_mean((i+1)/2));
                fprintf('\nROI Banana STD: %3.1f',roi.ban_sd((i+1)/2));
                fprintf('\nROI Coverage Percentage Compared To MRI volume: %2.1f%%', 100*roi.vol((i+1)/2)/mri.vol);
                fprintf('\nROI Banana Coverage Percentage Compared To MRI volume: %2.1f%%', 100*roi.ban_sum((i+1)/2)/mri.ban_sum);
                fprintf('\nROI Banana Coverage Ratio By Volume Compared To MRI volume: %2.0f:1 \n', roi.ban_sum((i+1)/2)*mri.vol/(mri.ban_sum*roi.vol((i+1)/2)));
            end
%             % Banana Stat - briks must be the same length to do stats.
%             stupid. even if were same dimension talairach transform makes
%             them different.
%             if roi.num > 1 &  % need to add in there that it must be the same dimensions
%                 roi_briks = cell2mat(roi_brik_reshape)';
%                 p = anova1(roi_briks);
%                 fprintf('\nANOVA between ROIs: p = %f',p);
%             end
        else 
            fclose(fid);
        end
        dig_brik = afni_dig_brik(choice,file,xyz,dig_brik,src,det,tai,opt_label,mri_dim,roi); % Creates the digitized brik with the appropriate colors and parameters
        write_brik(dig_brik,file,dir,exp,S,date);
%         choice.tai_banana_brik = input('Create TAI Banana Brik (Under Construction) (y/n) ','s');
%         if strcmpi(choice.tai_banana_brik,'y')==1
%             fprintf('\nUnder Construction\n');
%             % Can try using GCD method as long as it's not too large
%             % instead of recomputing from scratch
%             % mk_tai_banana_brik(choice,dig_brik,tai,xyz,opt_label,num_dum_pts,src,det,file);
%         end
	end
end

% Create and Write Thresh Scalp BRIK (Advanced mode only)
if choice.mode == 2
    choice.write_thresh = input('Create Scalp BRIK? (y/n) ','s');  
    if strcmpi(choice.write_thresh, 'y') == 1
        xyz_dig_brik = afni_dig_brik_ref_msk(xyz_mask);
        write_thresh(xyz_dig_brik,file,dir,exp,S,date,brik_suffix,thresh);
    end
end


% Create AOL file - an output ASCII file acpc space aligned data. 
% choice.aol = input('Create AOL File? (y/n) ','s');
% if strcmp(choice.aol,'y')==1
%     [xyz_aol] = mk_aol(xyz,dir,opt_label,exp,S,mri_dim);
% end
% 
% % AOL OUTPUTS: CAP STUDY ONLY. Outputs useful error information
% if strcmpi(exp,'cap') & strcmp(choice.aol,'y')==1
%     xyz_diff = mri_markers - xyz_aol(1:length(mri_markers),:);
%     for i=1:length(xyz_diff)
%         xyz_dist(i,1) = sqrt( (xyz_diff(i,1)^2) + (xyz_diff(i,2)^2) + (xyz_diff(i,3)^2) );
%     end
%     num_test_pts
%     avg_dist_tre =  mean(xyz_dist(3:length(xyz_dist)))
%     stdev_dist_tre = std(xyz_dist(3:length(xyz_dist)))
% end


% % Compute Map Error
% [xyz_aol] = mk_aol(xyz,dir,opt_label,exp,S,mri_dim)
% comp_map_error(file,mri_markers,xyz_aol)



% Transform the Data into a Talairach ASCII file
choice.tai_tfm = input('Talairach Transform Data(Create TOL File)? (y/n) ','s');
if strcmp(choice.tai_tfm,'y')==1
    mk_tol(choice,tai,xyz,opt_label,num_dum_pts,src,det,file);
end

% Remove Test, Ref, Param, Error from Subject and align2 folder
rm_align2_Sfiles(S_dbap,dir);


% Something is wrong w/ talairach BRIK - I give up
% % Write talairach BRIK 
% if strcmp(choice.tai_tfm,'y')==1
% 	tlrc_brik_choice = input('Write out talairach brik to view in Afni? (y/n) ','s');
% 	if (strcmpi(tlrc_brik_choice,'y')==1) & (strcmp(choice.tai_tfm,'y')==1)
%         tai.xyz = tai.xyz + repmat(tai.xyz(1,:), size(tai.xyz,1), 1);
%         xyz_tai = xyz_tai + repmat(tai.xyz(1,:), size(xyz_tai,1), 1);
%         tai.xyz = lpi2rai(tai.xyz);
%         xyz_tai = lpi2rai(xyz_tai); % loading back into RAI to view in afni
%         tlrc_brik_out_file = [brik_out_file '_tlrc.brk'];
%         dig_tlrc_brik = afni_dig_brik(xyz_tai,src,det,tai.fid,opt_label); 
%        	fid=fopen(tlrc_brik_out_file,'w');
% 		fwrite(fid,dig_tlrc_brik,'int8');
% 		fclose(fid);
% 		fprintf('Writing BRIK %s ... \n',tlrc_brik_out_file);
% 	end
% end

% save ([C:\Data\cap\ 'xyz' ] , 'xyz'); (saves one of the coordinate arrays
% to eventually make banana

% Closing Remarks
yodaman

%clear all

fprintf('\n\n  COREG COMPLETE\n\n');


